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Abstract 

The objective of this study is to propose a computational methodology that can effectively anchor the 
base flowfield of a four-engine clustered nozzle configuration. This computational methodology is based 
on a three-dimensional, viscous flow, pressure-based computational fluid dynamics (CFD) formulation. 
For efficient CFD calculation, a Prandtl-Meyer solution treatment is applied to the algebraic grid lines for 
initial plume expansion resolution. As the solution evolves, the computational grid is adapted to the 
pertinent flow gradients. The CFD model employs an upwind scheme in which second- and fourth-order 
central differencing schemes with artificial dissipation are used. The computed quantitative base flow 
properties such as the radial base pressure distributions, model centerline static pressure, Mach number 
and impact pressure variations, and base pressure characteristic curve agreed reasonably well with those 
of the measurement. 

Nomenclature 

A c area between nozzles in the exit plane 
A, vent area between nozzles 
C, turbulence modeling constant 

C 2 turbulence modeling constant 

G geometrical matrices 

h enthalpy 
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J Jacobian of coordinate transformation 

k turbulent kinetic energy 

M Mach number 

Pr turbulent kinetic energy production 

P pressure 

q represents l,u,v,w,h,k, and e 

S q source term for equation q 

t time 

U contravariant velocity 

u* wall-friction velocity, = (x w /p) 1/2 

u + non-dimensional velocity, = u/u” 

u,v,w mean velocities in x, y and z directions 
x,y,z physical coordinates 

y + non-dimensional distance, = y p u*p/p 

y p off-wall grid point distance to wall 

a conical nozzle half angle (=17.8°) 

p Prandtl-Meyer expansion angle 

y specific heat 

e turbulent kinetic energy dissipation rate 

p effective viscosity 

£, computational coordinates 

p density 

cr q turbulence modeling constants 

d> energy dissipation function 
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shear stress on the wall 


0 initial plume expansion angle 

grid adaption weighing parameter 

Subscripts 

a ambient or test cell 

b base 

be model centerline property on base 

e nozzle exit 

1 impact probe property 

o nozzle total property 

w solid wall 

Introduction 

Excessive base heating has been a problem for many launch vehicles 1 . For certain designs such as 
the direct dump of turbine exhaust inside and at the lip of the nozzle 2 , the potential burning of the turbine 
exhaust in the base region can be of great concern. Therefore, accurate prediction of the base environment 
at altitudes is very important during the vehicle design phase. Otherwise, the consequences could be 
undesirable. In the recent past, however, base environment of a launch vehicle has been predicted with 
large uncertainties using empirical methods, which either lead to out-of-database extrapolations, or overly 
conservative designs of the thermal protection system (TPS) and hence reduced payloads. The CFD 
method, which can be generically accurate when anchored, may provide a complementary prediction role, 
or be an optional design tool. 

In a previous study 3 , the turbulent base flowfield of a cold flow experimental investigation 4 for a four- 
engine clustered nozzle was numerically benchmarked. Parametric studies were performed on four 
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unadapted, relatively coarse algebraic grids (grid density varied from 34,030 to 113,202 points). 
Qualitative base flow features such as the reverse jet, wall jet, plume-to-plume recompression shock, and 
impingement have been captured. The physical nature of these flow features was in excellent agreement 
with that described in the experiment. Quantitative results such as the radial base pressure distribution, 
Mach number and static pressure variations along model centerline, were performed for a selected 
ambient-to-total-pressure ratio (P/P 0 ) of 39x1 0 J In addition, the base pressure characteristic curve was 
computed. These results agreed reasonably well with those of the measurement though relatively coarse 
grids were used. However, the trends of the model centerline Mach number and pressure distributions 
near the four-plume impingement point need to be improved, and the reversal ol the base pressure 
characteristic curve was not captured. Furthermore, the predictions for the radial base pressure and model 
centerline properties need to be broadened to other P a /P 0 ratios, and in general the prediction for the base 

pressure characteristic curve needs to be improved. 

Obviously, as pointed out in Ref. 3, the grid resolution played a dominating role in deciding the 
accuracy of the base flow solution. Higher grid density often resulted in better predictions. Also, when 
the grid lines that stemmed from the nozzle lip were specified at an angle corresponding to that ol a 
Prandtl-Meyer solution at P/P 0 = 39x10^, better predictions were obtained for the radial base pressure and 
the model centerline flow properties. In this study, further grid resolution studies were performed to 
demonstrate that the Prandtl-Meyer solution treatment for the initial plume expansion at different altitudes 
was highly efficient. In addition, as the solution evolved, flow gradient grid adaption was demonstrated 
to have greatly enhanced the efficiency and quality of the solutions, especially at higher altitudes where 
the plumes expanded to greater sizes and created stronger interactions. Pertinent base flow features such 
as the radial base pressure distributions, model centerline Mach number and pressure variations, and base 
pressure characteristic curve were computed and compared with the experiment, on a broadened database. 
Special base flow features such as the vent area choking and base shock were surveyed with the improved 
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solution. Knowing limited computational resources always prohibits unlimited increase of grid density; 
the proposed computational methodology will provide efficient and accurate base flow solutions for future 

launch vehicle TPS designs. 


Governing Equations 

The basic equations employed in this study to describe the base flowfield for a four-engine clustered 
nozzle are the three-dimensional, general-coordinate transport equations. A generalized torm of these 
equations written in curvilinear coordinates is given by 

(l/J)(3pq/3t) = 3[-pU,q + pC^Oq/d^)]/^, + (l/J)S q 

where q represents 1, u, v, w, h, k, and e, respectively. These are equations of continuity, x, y and z 
momentum, enthalpy, turbulent kinetic energy, and turbulent kinetic energy dissipation rate. The standard 
two-equation k-e turbulence model 5 closure is used to describe the turbulent flow. Turbulence modeling 
constants a q and source terms S q are given in Table 1. These turbulence modeling constants are widely 
used for nozzle flows 67 . 

The equation of state for an ideal gas is employed for the closure of the above system of equations. 
The characteristic of the governing equations changes from mixed parabolic-hyperbolic for subsonic flows, 
to mainly hyperbolic for supersonic flows. 

To solve the system of nonlinear partial differential equations, the methodology uses finite difference 
approximations to establish a system of linearized algebraic equations. An upwind scheme was employed 
to approximate the convective terms of the momentum, energy and continuity equations; the scheme is 
based on second and fourth order central differencing with artificial dissipation. The dissipation terms are 
constructed such that a fourth-order central and fourth-order damping scheme is activated in smooth 
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regions, and a second-order central and second-order damping scheme is used near shock waves. 

Viscous fluxes and source terms are discretized using second-order central difference approximation. 
A pressure-based predictor plus multi-corrector solution method is employed so that flow over a wide 
speed range can be analyzed. The basic idea of this pressure-based method is to perform corrections for 
the pressure and velocity fields by solving a pressure correction equation so that velocity/pressure coupling 
is enforced, based on the continuity constraint at the end of each iteration. Details of the present 
numerical methodology are given by ref. 6-7. 

Table 1 CT q and S q of the transport equations 
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Baseline Grid Generation 

A typical layout of an unadapted computational grid is shown in Fig. 1 . Due to the symmetrical 
nature of the flowfield, only 1/8 of this layout is generated and used for the actual calculation. The 
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boundary that is perpendicular to the center of the base is termed "model centerline". The two sides of 
the pie-shaped grid, as shown in Fig. 1, are the symmetry planes. For illustrative purposes, the symmetry 
plane that lies between the nozzles is named the "plume impingement symmetry plane", since the plume 
impingement line and the recompression shock will be on or attached to this surface, whereas the other 
one is termed the "nozzle symmetry plane", since it passes through the centerline of the nozzle. The 
centerline in which the two symmetry planes intersected is called the "model centerline . Two grid zones 
were created. The first zone started at the base and included the nozzle and the plume region. The second 
zone (the outer shell) is comprised of the ambient air, and a portion of the expanded plume. The baseline 
unadapted grids were generated using GENIE grid generator. The four nozzles, which are conical with 
a cylindrical external shell, are equally spaced on a circular base (heat shield) 4 , as shown in Fig. 2. The 
area ratio of the nozzles is 3.1 1 and the nozzle exit diameters are 2.67 inches. The base is located 2.0 
inches from the nozzle exit plane, giving a theoretical minimum vent area between nozzles of 
approximately 2.0 by 2.0 inches. The radial location of the theoretical minimum vent area, the four planes 
perpendicular to the base and between nozzles, is approximately 2.3334 inches from the centerline, which 
gives a vent area ratio (AJA C ) of approximately 0.96. This model is a larger scale model than the one 

used in Ref. 9. 


Boundary Conditions 

To start the calculation, an axisymmetric nozzle flow solution at the prescribed total pressure was 
carried out in a separate manner. A typical centerline exit Mach number for a total pressure of 60 psia 
was computed as 2.62. The converged flow solution was then mapped to a three-dimensional nozzle 
flowfield and the exit flow properties were specified as a fixed inlet boundary. The nozzle lip, nozzle 
outer wall and the base were specified as wall boundaries. The exit planes of zone 1 and zone 2, the outer 
surface (shell) of zone 2, and the inlet plane of zone 2 (flush with the base shield plane) were specified 
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as exit boundaries. In addition, a fixed (ambient) pressure was imposed on the inlet plane of zone 2, in 
order to obtain a unique solution for the corresponding altitude. Flow properties at the wall, symmetry 
plane, and exit boundary were extrapolated from those of the interior domain. 

A no-slip condition was imposed on the wall boundary and a tangency condition was applied at the 
symmetry plane. A modified wall function approach is employed to provide near-wall resolution which 
is less sensitive to the near-wall grid spacing. This is achieved by incorporating a complete velocity 
profile 10 . That is, 

u + = In [(y + + ll) 4 02 /(y + - 7.37y + + 83.3) 079 ] + 5.63 tan‘(0.12y + - 0.441) - 3.81 

This complete velocity profile provides a smooth transition between Logarithmic law-of-the-wall and linear 

viscous sublayer velocity distributions. 

Prandtl-Meyer Solution Treatment for Initial Plume Angle Resolution 
It has been shown 3 that the initial plume angle grid resolution is essential to the accurate prediction 
of base flow properties. In that study, the predicted base flow properties showed vast improvement even 
though a fixed initial plume angle (based on Prandtl-Meyer solution for P,/P 0 = 39x10^) was used. The 
natural extension of that work would be to construct an initial plume angle resolved algebraic grid for each 
pressure ratio according to the isentropic Prandtl-Meyer plume expansion theory. As shown in Fig. 3, the 
initial plume expansion angle can be expressed as 

0 = a + AP 

with 

(3 = [(y+1)/(y- 1)] 1/2 tan‘[(Y-l)/(Y+l) (M 2 -l)] 1/2 - tan '(M 2 -l) 1/2 
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where p is the Prandtl-Meyer expansion angle" through which a supersonic stream is turned to expand 
from M = 1 to M > 1. P, is based on the nozzle exit Mach number calculated from a simple one- 
dimensional calculation 12 . p a is based on the ambient-to-total-pressure ratio, which is equivalent to having 

a M a on the plume boundary. 


Solution-Adaptive Grid Generation 

A multi-zone, Self-Adaptive Grid Evolution (SAGE) code 13 , is used to refine the initial plume angle 
resolved algebraic computational grid. Its method is based on grid-point redistribution through local error 
minimization. The procedure is analogous to applying tension and torsion spring forces proportional to 
the local flow gradient at every point and finding the equilibrium position of the resulting system of grid 
points. Since Mach number contour is closely associated with the plume boundary layer, whereas the 
pressure gradient follows the recompression shock, these two flowfield variables were used as pertinent 
grid adaption parameters. The adaptive function is a combination of both and can be expressed as 

C M aM/ac + C P 3p/ 3$ = + ^ p s 

Fig. 4 shows slices of four typical computational grids. Each slice is a portion of the nozzle symmetry 
plane and is bounded by the nozzle centerline and the model centerline. Grid A is an algebraic grid 
treated with Prandtl-Meyer solution for the initial plume angle resolution. Grid B is the result of Grid 
A adapted solely to a pressure solution. The clustered grid lines clearly exhibit the plume-to-plume 
recompression shock, although the shock on the nozzle symmetry plane is not as strong as that in between 
two nozzles, or the "minimum vent area" plane. Grid C is the outcome of grid A adapted entirely to a 
Mach number solution. The packed grid lines follow the plume boundary and the initial plume expansion 
angle resolved algebraic grid lines that stem from the nozzle lip. Notice the adaption was applied several 


9 


grid lines above the nozzle lip so as to maintain the initial expansion resolution. Grid D is the adaption 
of Grid A where 50% Mach number gradient and 50% pressure gradient were used as the adaptive 
function. The grid line clustering follows both the plume boundary layer and the recompression shock. 

The computations were performed on a NASA/MSFC CRAY-YMP. The computational time lor a 
typical calculation was estimated as 1.0x10^ CPU seconds per grid per step. Approximately 3000 to 4000 
iterations were required for a 119,016 grid points solution to reach approximate convergence and an 
additional 2000 iterations were needed for a higher grid density (e.g., 168,399 grid points) solution to 
converge when the initial flowfields were started afresh. The storage requirement of the CFD model is 

40 words per grid point. 


Results and Discussion 

Static Pressure. Mach Number, and Impact Pressure Variations along Mo del Centerline 

Static pressure, Mach number, and impact pressure comparisons along model centerline, assess the 
accuracy of the model prediction for the strength of the reverse jet. The accuracy of the model prediction, 
however, depends on proper computational grid distributions. Fig. 5 shows a comparison of variations 
along model centerline for P a /P 0 = 39x1 0 4 . The effects of the Prandtl-Meyer solution treatment for initial 
plume angle resolution and solution-adaptive gridding were obvious: the solution with the highest grid 
density (245,493 grid points) in which initial plume expansion was not resolved, although it employed 
twice the number of points in the initial plume boundary layer (in comparison to the 1 19,016 points setup) 
and had the grid refined according to the pressure gradient, produced the worst comparison; whereas the 
solutions using a 45. 4-degree initial plume angle resolution, including one solution that ran without any 
grid adaption, computed significantly better agreements with less than half of the grid points. Among 
those three solutions, the one using pressure-solution gradient adaption produced best centerline property 
comparisons with those of the experiment. The one using Mach number-solution gradient adaption did 
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not compare as well. This is because the model centerline properties are enclosed by the recompression 
shock where the solution is best predicted by grid following pressure gradient. This example clearly 
demonstrates the validity and efficiency of the proposed computational methodology. 

Fig. 6 show the comparisons of static pressure, Mach number and impact pressure variations along 
model centerline for three ambient-to-total-pressure ratios. In subsonic flow region, the impact pressure 
was reduced through 


p. = P (1 + (y - l)/2 M 2 ) T ' (rl) 

In supersonic flow region, assuming the pitot tube 4 is immersed behind a shock wave, the Rayleigh pitot 
formula 11 was used. 

P. = p [(y+ l)M 2 /2]^ 1) {(Yfl)/[2YM J -(Y- 1 )] } Wrl> 

Although extensive grid parametric studies have been performed, for the purpose of clarity, only selected 
comparisons were shown. It can be seen that different initial plume expansion angles were applied for 
different ambient-to-total-pressure ratios. At P a /P 0 = 20x10"*, where maximum model centerline base 
pressure and peak Mach number were measured, 168,399 grid points were required for additional 
resolution. In general, The predictions agreed reasonably well with those of the experiment. The impact 
pressure decreases from the plume impingement point (approximately at Z = 4.5 inches) to the reverse jet 
recompression location due to radial flow. Downstream of the recompression it remains constant due to 
the prevailing subsonic flow. In general, the peak Mach number increases and the valley static pressure 
decreases as the pressure ratio (P/PJ drops, and the position of the peak Mach number moves toward the 
base as does the valley of the static pressure. The strength of the reverse jet also increases as the pressure 
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drops. At P/P 0 = 20xl(T\ a base shock is formed between 0.5 to 1.0 inches from the base center. Due 
to high viscous dissipation and in general a weak reverse jet (in comparison to under-expanded supersonic 
nozzle jet direct impingement on a perpendicular surface), the base shock is a smeared shock 4 , as 
evidenced by the moderate increase of static pressure over a finite distance. Fig. 7 shows the computed 
iso-pressure surfaces at P/P 0 = 20x1 O' 4 . The computed plume-to-plume recompression shock (iso-value 
= 40 lb/ft 2 ) resembles closely that of a S-IV four engine stage exhaust plume Schlieren photograph 14 . The 
base shock (iso- value = 18 lb/ft 2 ) is situated above the heat shield, in between nozzles. 

Radial Base Pressure Distribution 

Radial base pressure data were taken at the base along the plume impingement symmetry plane, hence 
the comparisons benchmark the model for the predictions of the reverse jet at base center and the wall jet 
in the vent area. It can be seen from Fig. 8, the computed radial base pressures agreed reasonably well 
with those of the experiment. For all three ambient-to-total-pressure ratios, the reverse jets had formed 
and the peak pressures occurred at the base center, whereas the radial base pressure decreased as the 
distance from the center of heat shield increased. The radial base pressure eventually dropped to that of 
the cell pressure, which is physically correct. 

Base Pressure Characteristic Curve 

The center base pressure variation with ambient pressure (altitude) has become known as the 
characteristic curve 4 . Representing the location of the severest environment on the base, it is one of the 
important parameters in designing the thermal protection system for the launch vehicles. Fig. 9 shows 
a comparison between the measurements and the predictions. Matz and Goethert data from which an 
identical AJA C ratio of 0.96, albeit a mere 0.80 inches distance between the base and nozzle exit plane, 
were selected and plotted along with Brewer data 4 for background comparison. Both experiments were 
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conducted over the characteristic range of Pbc/P a from near 1.0 to near 4.0, corresponding to altitudes 
ranging from 22,800 feet to 122,500 feet. At ?J? 0 = 100x10" or PJP a = 1, the four exhaust plumes do 
not interact much with each other, except some aspiration exists. As altitude increases or as pressure ratio 
drops, the plumes start to interact more and the aspiration decreases; In the mean time, base pressure 
decreases whereas a reverse jet and subsetjuently the wall jets take shape. The predicted base pressure 
characteristic curve agreed very well with those of the experiment 4 . The condition of Pbc/P, = 2 indicates 
a choking condition for the wall jet if the system was an enclosed isentropic convergent-divergent nozzle, 
which obviously does not apply since we are dealing with a complicated three-dimensional turning wall 
jet in a vent area where the open-top can only be closed by expanding sonic plumes. The fact that Mach 
number gradient adaption was applied mostly for P/P 0 > 26x10" whereas pressure gradient adaption was 
used at lower pressure ratios indicating that characteristic base pressure is dominated by plume boundary 
layer resolution when the plumes are further apart; as the plumes close, the recompression shock becomes 
more important, and better base pressure was resolved by applying P 5 adaption. In fact, the plumes 
closed completely at approximately P/P„ = 23x10" where the minimum base pressure occurs. After that 
enclosure, P* raises to its maximum at P/P 0 = 20x10" where the wall jet boundary layer grows and 
accelerates to sonic velocity in the vicinity of the vent area at P,*/P, = 4. P^ then decreases as the ambient 
pressure continues to drop until the vent areas are completely choked. Further reduction in ambient 
pressure would not change the base pressure significantly after the total choking, as indicated by the 
leveling-off of the predicted P* from P/P 0 = 10x10" to 1x10", albeit the experiments 49 stopped at P a /P 0 
= 15x10" and 10x10", respectively, due to hardware limitations. 

Fig. 10 shows the computed sonic surfaces at P./P 0 = 1x10". There is an approximate ellipsoid-shaped 
sonic surface in the middle of the domain, manifesting the acceleration of the reverse jet. The complex- 
shaped vent area sonic surface, not resembling the theoretical minimum vent area plane (Fig. 2), sealed 
all the flow path with the plume sonic surface and created the total choking. 
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Conclusion 


A computational methodology has been developed to benchmark the base flowfield of a four-engine 
clustered nozzle configuration. It is based on: a three-dimensional, viscous flow, pressure-based CFD 
formulation, a Prandtl-Meyer solution treatment for the algebraic grid which is proved to be 
computationally efficient for the initial plume expansion resolution, and the computational grid which is 
subsequently refined according to pertinent base flow physics. The predicted physical flow features such 
as the reverse jet, wall jet, recompression shocks due to plume-plume and reverse jet-base impingement, 
plume enclosure, and vent area choking are in reasonable agreement with those described in the 
experiment. The predicted quantitative results such as the radial base flow distribution, static pressure, 
Mach number and impact pressure variations along model centerline, and the base pressure characteristic 
curve also agreed well with those of the measurement. This methodology not only provides insight into 
the multiple engine base flow physics, but also will be useful in the design and analysis of TPS for launch 
vehicles. 
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Fig. 3 Prandtl-Meyer expansion around the nozzle lip. 



Fig. 4 Slices of four different computational grids. 
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Fig. 5 Comparison of Mach number and static pressure variations along model centerline 
for P JP D = 39x1 0 4 
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Fig. 6 Comparisons of static pressure, Mach number, and impact pressure variations 


along model centerline. 
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Fig. 8 Comparisons of radial base 







Fig. 9 Base Pressure Characteristic Curve. 
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